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ABSTRACT 

Minimum variance estimation requires that the statistics 
of random observation errors be modeled properly. If 
measurements are derived through the use of atomic 
frequency standards, then one source of error affecting 
the observable is random fluctuation in frequency. This 
is the case, for example, with range and integrated 
Doppler measurements from satellites of the Global 
Positioning System used for precise geodetic point 
positioning and baseline determination for geodynamic 
applications. In this paper an analytic method is 

presented which approximates the statistics of this 
random process. The procedure starts with a model of 
the Allan variance for a particular oscillator and 
develops the statistics of range and integrated Doppler 
measurements. A series of five first order Markov 

processes is used to approximate the power spectral 
density obtained from the Allan variance. Range and 
Doppler error statistics are obtained from the integra- 
tion of the corresponding autocorrelation function. 
Statistics for residuals to polynomial clock models are 
then obtained by linear transformation. Examples are 
given for rubidium and cesium clocks. 


ATOMIC CLOCK ERRORS AND FREQUENCY STABILITY 

A clock is any device which counts the cycles of a periodic phenomenon 
and among the most stable clocks in use are the atomic clocks which 
form the basis for atomic time scales such as International Atomic 
Time (TAI). Atomic time is used primarily as a measure of time inter- 
val and is based on the electromagnetic oscillations produced by quan- 
tum transitions within the atom. The precise definition of stability 
is found in Blair (1974). Basically it is a measure, usually given 
statistically, of the random fluctuations in frequency which can occur 
in a clock's oscillator over specified periods of time. For a given 
time interval a particular oscillator is considered best if the ex- 
pected level of frequency fluctuation is a minimum in terms of the 
Allan variance defined below. 

Equation (1) is the model used to describe the types of error present 
in atomic time scales 


551 


( 1 ) 


Ti(t) = i D^(t - t^)2 + R^(t - t^) + T^(t^) + x(t) 

The deterministic errors consist of bias, drift, and ageing terms 
modeled as a quadratic polynomial in time. The ageing term is less 
observably for clocks whose long-term stability is good such as cesium. 
The term x(t) in equation (1) represents the random time error due to 
the integration of random fluctuations in frequency: 


x(t) = j- S f(x)di = / y(l)dT 
It t 


( 2 ) 


The magnitude of this term depends on the stability of the clock and on 
the interval of time which has passed since the scale was reset or 
calibrated. 


Hellwig (1977) points out that "the characterization of the stability 
of a frequency standard is usually the most important information to 
the user especially to those interested in scientific measurements and 
in the evaluation and intercomparison of the most advanced devices 
(clocks)." Since the frequency stability of a standard depends on a 
variety of physical and electronic influences both internal and exter- 
nal to the standard, measurement and characterization of frequency 
stability are always given subject to constraints on environmental and 
operating conditions. In addition frequency stability depends on the 
exact measurement procedure used to determine stability. 


Frequency stability characterization is done in both the frequency and 
time domain. In the time domain a frequently used measure of stability 
is the Allan variance or its square root. In the frequency domain it 
is the power spectral density. 


The Allan variance as a time domain measure of frequency stability is 
found especially useful in practice since it is obtainable directly 
from experimental measurements and contains all information on the 
second moments of the statistical distribution of fractional frequency 
error. The Allan variance is defined as follows: let y. 


^k’ ^k+l’ 


'k+2’ 


^ 1 ’ ^ 2 ’ • • 

be observed fractional f requency^errors separ- 


ated by a repetition interval of T seconds. For each integer N greater 
than or equal to 2, calculate y^ from 


^m = N 


(m+l)N - 
I 

k = mN 


m = 0, 1, 2, 


M. 


(3) 
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The Allan vari- 


This is an average over N consecutive values of y . 

2 - * 
ance, a^(N), is then obtained from the averages by 




M - 1 
1 

m = 0 


<Vi - 


(4) 


An examination of this equation reveals that the Allan variance for a 

particular sampling interval NT is the average two-sample variance of 

the y (N) . 

•'m 

For frequency standards the square root of the Allan variance is usual- 
ly given in graphical form on a log-log scale. For individual classes 
of frequency standards models for the Allan variance are used which 
portray general frequency stability characteristics. Hellwig (1975) 
gives examples of such models for many oscillator types. Figure 1 
shows the typical form. In this form, a (x) is the square root of the 
Allan variance for the sample interval t.^ The quantity is called 
the flicker floor and Tj, X-, are the break points of the plot. The 
constants associated with this figure are usually specified for each 
type of frequency standard. A comparison of such information can 
facilitate the selection of a frequency standard for a specific appli- 
cation. 


The stability characteristics shown in the three regions of Figure 1 
are typically present in many Allan variance plots of specified oscil- 
lator performance. The first part, region I, reflects the fundamental 
noise properties of the standard. This behavior continues with in- 
creased sampling time until a floor is reached corresponding to region 
II. After the performance deteriorates with increased sampling 
time. Hellwig (1977) outlines the error sources corresponding to each 
portion of the graph. The magnitude and slope of each segment will 
depend on the particular category of standard. 


An alternative procedure for specifying the stability of a frequency 
standard, in the frequency domain, is the use of the power spectral 
density (PSD) of instantaneous fractional frequency fluctuations y(t) . 
Allan et al. (1974) have given a useful model to represent the PSD for 
various categories of frequency standards. This model is in the form 
of a power law spectral density having the form 


s 

yy 


(u)) = 




0 ^ u) ^ 




U) > 


(5) 


where a takes on the integer powers between -2 and 2 inclusive depend- 
ing on how the interval (0,U)j^) is to be divided into subintervals, one 
for each a to be used. The quantity h is a scaling constant, and the 
PSD is assumed to be negligible beyond the frequency range (0,u».). 
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LOG iOy) 


I. Qy - t ’ or T ’/2 


II. Qy = CONSTANT 

III. Oy-T° 0 < a < 1 


TIME (SEC) 


LOG (t) 


Fig. 1-General frequency stability characteristics 

Barnes et al (1971) and Meditch (1975) give the transformations between 
the time domain measures of frequency stability in the form of the 
Allan variance and the power law spectral densities. Table 1 taken 
from Meditch gives these conversions for three types of fractional fre- 
quency error sources. 

Table 1-Allan variance and power spectral 
density for common error sources 



ERROR SOURCE 


WHITE NOISE 


FLICKER NOISE 


INTEGRAL OF 
WHITE NOISE 
(RANDOM WALK) 


ALLAN VARIANCE 

(T) 


TWO SIDED SPECTRAL DENSITY 

Syy (m) 
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RANGE AND DOPPLER OBSERVATION ERRORS DUE TO RANDOM ATOMIC CLOCK ERROR 


As previously discussed, an atomic clock's time scale can be expected 
to differ from ideal time due to both deterministic and random errors. 
The random component is due to integration of fractional frequency 
errors. A range observation determined from radio signals broadcast by 
a satellite is subject to the random errors of the frequency standards 
in both the satellite and the tracking receivers. The effective range 
error at time t due to the timing error in one of the time scale is 

6R.(t) = cT.(t) (6) 

with the random component being the random walk 

t 

n^(t) = c / y(t)dT (7) 

^s 

where c is the velocity of light. The random component is due to the 

accumulated effect of fractional frequency error since the clock's 

start or reset at t . 

s 

The random error r|-(f) is correlated in time. Consider two measure- 
ments of range R(t.) and based on the use of the oscillator in 

the satellite, and ilssume momentarily that the receiver's oscillator is 
free from random error. The covariance between these measured ranges 
due to correlated fractional frequency error in the satellite oscil- 
lator is 


E[R(t )R(t )] = E[n(t )n(t^)] 

V J 

2^1 ^k 

= c^E[J J y(T)dT ; y(T')dT^] 

t t 

s s 

= c ; ^ / E[y(T)y(T')]dtdt' 


where d) (t - T^) 

yy 

quency error y(t) 



J 4» (t - t')dxdT ' 
^ yy 


( 8 ) 


s s 

is the autocorrelation function for fractional fre- 
defined by 
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♦yy(x - x') = E[y(t)y(t')] 

00 00 

= / J yy'f(y,y',t,t')dy dy'. (9) 

-00 -00 


The function f(y,y is the joint probability density function for 

fractional frequency error. Here it is assumed that y(t) is a mean 
zero stationary random process. The function “ X could be 

obtained by the inverse Fourier transform of the given power > spectral 
density S (u)) : 

yy 




00 

/ S (u))e^***^du) 

-i yy 


( 10 ) 


where 


t = X - x'. 


An alternate procedure for obtaining the autocorrelation function 
(b (t) from the Allan variance is given below. 

yy 

The variance of a range observation is obtained from equation (8) by 
setting tj equal to t^^: 




(j) (x - x')dxdx'. 
yy 


( 11 ) 


The presence of random frequency error in the receiver oscillator 
introduces additional, but similar, terms into equations (8) and (11) 
which must be considered when assessing the range uncertainty due to 
all random clock errors effecting the measurement. 


For integrated Doppler or range difference observations the random 
measurement error associated with system clocks is the integral of 
fractional frequency error over the Doppler integration interval. The 
random error in range difference due to one oscillator is 

Hij = n(tj) - n(t^) 

t. 

= c / ^ y(x)dx. (12) 

t. 
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Notice in equation (12) that the random error n. . is a function of t., 

ij i’ 

t., and y(t). The error does not depend on t . Range difference 
J s ° 

measnrements have the following correlation for each oscillator 


E[Mi_AR^] = E[n,.n^^] 

t_. t. 


= J ^ ^ - x')dTdx 

t. 

1 k 


with the variance 


^AR. . 
xj 


, t. t. 

= c / J J A (X - x')dxdx'. 


t. t. 
1 1 


yy 


(13) 


(14) 


Observe that the random range difference errors, whose statistics are 
given by equations (13) and (14), are stationary; however, random range 
errors, whose statistics are given by equations (8) and (11), are not, 
A stationary random process is one whose statistics are invariant in 
time . 


For the oscillator performance specifications shown in Figure 2 exam- 
ples of the contribution to the range error are given for both oscilla- 
tors in Figures 3 and 4 over a five-day span. The clocks are assumed 
to be perfect initially. Also included is the standard error for the 
random walk n(b) obtained using equation (11). The procedure used in 
simulating the random range error is discussed in Meditch (1975). 



Fig. 2-Allan variance for satellite and receiver oscillators 
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RANGE QUANTITIES IMETERSI 



TIME (DAYS) 

Fig. 3-Standard error and random range error 
based on receiver cesium specifications 




RANGE AND DOPPLER OBSERVATION ERROR STATISTICS 


Fractional Frequency Autocorrelation from the Allan Variance 

The equations giving the second order statistics of random range and 
integrated Doppler observation errors due to random fractional fre- 
quency errors were presented in the last section. Those equations 
require that the fractional frequency autocorrelation function be 
known. In this section discussion of a procedure for obtaining an 
analytic approximation to this function from the Allan variance is 
given. This method yields a simple analytic autocorrelation function 
and avoids numerical difficulties that may arise when the inverse 
Fourier transform of the power spectral density is evaluated. 

The Allan variance models shown in Figure 2 for the satellite rubidium 
and receiver cesium oscillators are a function of the sampling time T 
having the form 


!!o 

T 


< T < T, 



Using the transformations in Table 1 the power spectral density for 
fractional frequency may be developed from equation (15): 


S 

yy 


(u)) 



(16) 


The square roots of the power spectral densities corresponding to the 
Allan variance specifications of Figure 2 are given in Figure 5. The 
constants associated with the two functions and the formulas for 
computing the constants associated with the power spectral density 
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function based on the Allan variance are given in Table 2. These for- 
mulas are developed from the transformations of Table 1. 

The autocorrelation function <)>yy(t) can be obtained from the power 
spectral density using equation (10) 

^ Syy(u))e^‘“^du). 

However, as a result of transforming the band limited white noise 
portion of the spectrum, this form for the autocorrelation function has 
an oscillatory behavior for small t. This is an artificiality of the 
model . 



Fig. 5-Square root of PSD 


An alternate approach for obtaining an autocorrelation function is to 
approximate the power spectral density model with a smooth function 
whose autocorrelation is expressible in simple analytic form. The 
first step in this development is to approximate the flicker noise 
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Table 2-Oscillator parameters 


QUANTITY 

UNITS 

FORMULA 

SATELLITE CLOCKS 
(RUBIDIUM) 

RECEIVER CLOCK 
(CESIUM) 


sec 


I.OOxlo’ 

1.00x10® 


sec 


IMxlo® 

1.00x10® 

^3 

sec 


1 .(»xio« 

1.00x10' 

“0 

sec"^ 

yiiiTj 

1.73x10‘ 

1.73x10“' 

0), 

sec"’ 

6 ln2/(nT2) 

1.32x10 ® 

1.32x10“' 

Ci>2 

sec"’ 

n/(2T,»n2l 

2.27x10 ’ 

2.27x10® 

Of 



6.00x10^” 

3.00x10“’* 

No 

sec 

t,o,2 

3.60x10 “ 

9.00x10"“ 

N, 


noj^/2 In 2 

8.16x10 “ 

2.04x10“" 

Nj 

sec"’ 

3of^h2 

1.08x10 ” 

2.70x10““ 

N 3 

sec 


3.60x10 

9.00x10““ 

a 



2.36x10® 

1 61x10® 


sec"’ 


2.03x10 ® 

1 67x10 ® 


segment of the spectrum by a series of cascading functions whose values 
alternate between being constant and being inversely proportional to 
the square of the frequency. This type of procedure is described by 
Meditch (1975) in constructing a linear system which simulates flicker 
noise using a white noise input. Figure 6 shows the transfer function 
for flicker noise. A three stage cascading transfer function is super- 
imposed consisting of the functions F , Fg, and F^ wMch are defined in 
Table 3. These functions are defined to have the required properties 
and give a continuous although not smooth approximation to the flicker 
noise power spectral density. 

The constants of this approximation are now derived over frequency 
intervals as given in Meditch (1975). The general form of the function 


Fa(uj) = (17) 

a 

between the frequencies u) and oiu . At frequency uj , defined in 

Table 2, the function F. takes on the value 

A 


561 


SQUARE ROOT OF PSD 



Fig. 6-Three stage transfer function approximation 
of flicker noise spectrum 

Table 3-Definition of three stage transfer function approximation 

FUNCTION INTERVAL DEFINITION (PSDI 



^ U> ^ (Uf 

N,/u, 


CO, < Cij ^ aw. 

N*/o,2 


aw, < w ^ w^ 



w^ ^ w ^ a^w. 

N^lo^w,^ 


o^w, < w < a^w. 

Ng/w^ 


a^w, < w < wg 

Nb/o‘c-.» 


Wg < w ^ a^w. 

NgloOu.J' 

fc 

a^w < w < a*w. 

Nc/u,2 

WHERE 

a®w < w ^ ui2 

N„ 

1 


- aw^N^ 



Nb = o^u-iN, 

w, = W,\/T 


Nc = o5w,N, 

n = 3 
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( 18 ) 


U) 

a 


N 

lU 


1 

1 


since the flicker noise power spectral density has the same function 
value at frequency uu^ . Solving equation (18) gives 


N u)^ 


(19) 


A similar analysis gives the constant N . The function F has the form 

D D 


N. 


B 


FbW = 2 

UJ 


( 20 ) 


At frequency o lu^, Fg has the function value 


2 

''b‘“ “a’ = -Tl 


N, 


2 2 
a lu a u) 
a a 


( 21 ) 


since at a u) the function F„ has the same value as function F. at 
3 15 A 

frequency Cdu (see Figure 6). Solving equation (21) and using equa- 
tion (19) ^ 


= a\ N,. 

B A 11 

For the function F^, 


( 22 ) 


Fc(u.) = -f 

UJ 


(23) 


its function value at frequency a u)^ equals the value of Fg at fre- 
quency a^u) giving 


‘“a^ 8 2 6 2 

a u) a u) 
a a 


(24) 


Using equation (22) gives the solution 

"c = “'"b = “S"!- 


(25) 
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Numerical values for a and u) are given in Table 2. The power spectral 
density consisting of the tSree cascading functions and the remainder 
of the original function will be denoted as the second power spectral 
density model for each oscillator. 

The next step in the development of a simple analytic autocorrela- 
tion function is to approximate various segments of this second model 
with a first order Markov process power spectral density function, a 
function of the form 


S(iu) = 


2 

2a 3 

2 ^ q2 

lu + 3 


(26) 


where 3 is the inverse of the correlation time (see Gelb (1974)). The 
autocorrelation function for a first order Markov process is 


4»(t) = a^e"^ ^ . (27) 

Notice in equation (26) that the power spectral density decreases as 
the inverse of the square of the frequency. This is the type of func- 
tional behavior seen in the interior of the cascading functions 
through F„. It is also the behavior of the original power spectral 
density in the interval (o»q, lu^). In addition the power spectral 
density of the Markov process remains virtually flat until the fre- 
quency reaches a point at which the function decreases rapidly. These 
properties make this function an excellent choice for approximating the 
second power spectral density model piecewise. 

The second model is then divided into five segments defined in Table 4. 

The high frequency cut off shown as 10 ^ in Figure 5, will be in- 
creased so that the band limited white noise component of the power 
spectral density may be approximated better by the first order Markov 
power spectral density. 


Table 4-Division of second PSD model for Markov process approximation 


NOTATION 

•l 

*2 

*3 


INTERVAL 

(0, cu,] 

(o)i, acoj 
law,, a^w,] 
[o^w,, a®w,] 
[a^w,. w^J 
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The approximation consists then of fitting a function in the form of 
equation (26) to each subdivision of the second model S' (uj) given in 

Table A. There are two parameters a and p to be determined for each 
segment giving a total of ten parameters. 


The procedure which was adopted was an asymptotic approximation whereby 
two constraints were imposed on the Markov power spectral density 
function giving a and ^ directly. This procedure was implemented 
because of simplicity and because the results compared favorably with a 
least squares approach. The asymptotic approach develops an approxima- 
tion on the interval 1., 



using the following constraints: 


(i) at zero frequency the approximating Markov power spectral 
density equals the second model at frequency 

S.(o) = S^(u)^) (28) 

(ii) in the limit as u) increases the value of the function S . (u)) 
converges to the following function 


lim S (u)) = 

U) -► e»J U)^ 


(29) 


and at u). this limiting value is set equal to the value of S' (u») : 
X yy 



s' (U)-). 
yy A 


(30) 


Equations (28) and (30) are a system of two equations in two unknowns. 
Their solution yields the parameters a. and p. for the approximating 

Markov power spectral density function S . (lu) . The nature of the second 
constraint, equation (30), is to force the function Sj (u)) to asymptoti- 
cally approach S^(u>) at The first constraint is necessary to ap- 
proximate the white noise or flat component of S' (w) at the beginning 
of each subinterval. ^ 


Finally a comment concerning the approximation in the last subdivision 
is necessary. In order to obtain a good approximation to S^(u>) in 

that interval it is necessary to choose u), large enough to allow the 
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flat portion of the Markov process spectral density to fit the white 
noise component which dominates this interval (see Figure 5). Choosing 
three or four orders of magnitude larger than 0.1 and two or 

three orders of magnitude smaller than enables a good approximation 
to be made but adds power at these higher frequencies. The result is 
an autocorrelation function which tends to a delta function as u)^ 
goes to infinity and whose variance increases as in is chosen larger 
(see Figure 7). However, this will have negligible effect on range 
and range difference statistics. 


The smooth fractional frequency autocorrelation function 

given by the inverse Fourier transform of the five Markov process power 
spectral densities S . (u>) . The result of each transformation is an ana- 
lytic function whose'^form is given by equation (27). The final result 
is the sum of these functions 



(t) 



2 

a.e 

J 




(31) 


For range and integrated Doppler observations the statistical contribu- 
tion due to random oscillator error is obtained using equation (31) in 
equation (8) through (14). 

Figures 8 and 9 show the original transfer functions and the asymptotic 
approximations. The parameters obtained using this approximation 
procedure are given in Table 5. 

Observation Error Statistics Based on Markov Process Approximations 

The first order Markov autocorrelation function, equation (31), and 
equations (8) through (14) give the second order statistics for random 
range and integrated Doppler observation errors due to each oscillator 
used in the measurement process. These integrals may be evaluated 
giving analytical expressions for the variance and covariance of range 
and Doppler observations. 

Let R(t^) and R(tj^) be range observations subject to one random clock 
error only. The covariance between the observations is given by equa- 
tion (8). Using the first order Markov approximations, the integration 
of equation (8) gives the covariance as 
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Fig. 9-Receiver cesium transfer function and sum 
of asymptotic approximations 

Table 5-Fractional frequency autocorrelation function 
parameters for Markov process approximations 


ASSYMPTOTIC LEAST SQUARES 


OSCILLATOR TYPE 

INTERVAL 

(ALPHA)2 

BETA 

(ALPHA)* 

BETA 

RUBIDIUM (SPEC! 


3.1177x10“” 

1.732x10"* 

3.3719x10"** 

1.681x10* 


>2 

6.2625x10"“ 

2.032x10"' 

7.6238x10"“ 

2.256x10"' 


<3 

6.2625x10"“ 

1.128x10"* 

7.6246x10"“ 

1252x10 * 


<4 

62625x10"“ 

6262x10 * 

7.6246x10"“ 

6.947x10* 



1.8000X10"’* 

1.000x10** 

1.9343x10"’* 

9.631x10** 

CESIUM (SPEC) 


7.7942x10"” 

1.732x10"* 





12922x10"” 

1.677x10 * 




■3 

12922x10"” 

4221x10"' 




•4 

12922x10"” 

1.113x10"' 





4.5000x10"“ 

1.000x10** 




- 1.0x10* S - Nq/I OxIO* 
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E[R(t.)R(tj^)] 


E(n(t^)n(tj^)] 

5 5 

J = 1 L 


- S’ * r ( ' 

J 




+ e J * ® - e 


-Sj‘s - S’ . 


( 32 ) 


for tj^ greater than t., where t is the start or reset time of the 
clock. The variance of the random range error is obtained by setting 
tj^ equal to t^ in equation (32) 


E[R(t.)R(t.)] = E[n(t.)n(t.)] 





(33) 


The range error r|(b) resulting from the integration of fractional 
frequency error y(t) is a statistically nonstationary process. An 
examination of equations (32) and (33) reveals terms which are func- 
tions of t., or tj^, minus t . Thus, for instance, the variance in- 
creases wiA time. This is vllustrated in Figure 10 for the rubidium 
clock. The standard error of a range measurement based on the use of 
this clock is given for 20 range observations spaced at 15-minute 
intervals starting five minutes, one hour, and five hours after the 
start of the clock. The increase in variance is almost linear. An 
examination of the autocorrelation function shows that this function, 
dominately flat, is similar to a random bias having a constant auto- 
correlation and whose integral is a random ramp which increases exactly 
linearly. Hence a linear growth in variance is expected as seen in 
Figure 10. The correlation coefficients between the first and the 
i'th range observation in each of these sequences are given in Fig- 
ure 11. As the starting time of the sequence increases from t , so 
does the correlation among the random errors. This again is expected, 
since the variance increases with time and the errors are correlated. 

Figure 12 gives the autocorrelation function for the cesium clock based 
on the Markov process approximation and Figures 13 and 14 give the 
standard error and correlations of range errors based on this clock. A 
comparison of Figures 10 and 13 reveals the greater stability of the 
cesium clock. After ten hours of operation the standard error of the 
cesium clock output is approximately 3.5 nanoseconds compared to 63 
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s '2 15 ig 21 

NUMBER OF 16 MINUTE RANGES 


24 


Fig. 11-Correlation coefficients between range 1 
and range i (rubidium clock) 
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Fig. 14-Correlation coefficients between 
range 1 and range i (cesium clock) 

nanoseconds for the rubidium standard. In addition, the correlations 
among the cesium clock errors decrease more rapidly than the rubidiiun 
clock errors. Considering both random clock error sources the total 
variance and correlation of range observations Rj^(t.) and Rj^(t.) meas- 
ured by receiver k are given by the equations ^ 

E[R^(ti)Rj^(ti)] = E[n^(t.)n 3 (tp] + E[n^(t.)n^(tp] (34) 


E[\(ti)R^(tj)] = E[ri3(t.)ng(t.)] + E[rij^(tpnk(tj)] 


where the variances and correlations of the random error r| are given by 
equations (32) and (33) . The subscript "s" refers to the satellite 
rubidium clock. 


For simultaneous observations of range by two receivers the covariance 
of the observations Ejj(b^) and R^^(tj) is given by 

E[R^(t.)R^(t.)] = E[ng(t.)n^(tj)]. (36) 

In the above equations the random errors n have zero mean which is a 
consequence of fractional frequency error being zero mean. 

Let AR(t ) be an integrated Doppler or range difference measurement 

over the interval (t., t ) and AR(t„) a similar measurement from the 

in X, 
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same receiver over the interval (tj^, The covariance of the obser- 

vations is 


E[AR(t^), AR(t^)i = Ein(t^) - n(t.), n(t^) - 

= E[n(t^)n(tj^)i - E[n(t^)n(tj^)] - E[n(t.)n(tj] 
+ E[n(t.)n(t. )] 


= c 


i' - 'k 




-p (tj - t„) - t„) 

e - e 


-p (tj - -Pj(^ - 1.) 

- e '* + e 


(37) 


The variance of a range difference observation is given by 


5 2a^ ^ 


E[AR(t )AR(t )] = c" I 

n n ■ _ 1 p- 

J - 1 J 


1 


J 




. (38) 


Equations (37) and (38) are independent of the clock epoch t^. The 

statistics of the range difference error depend only on the Doppler 
integration interval or the time difference between observations. Thus 
the random range difference error is stationary. Expressions analogous 
to equations (34) through (36) express the complete statistics of range 
difference observation errors for individual or simultaneous observa- 
tions due to clock error. 


STATISTICS OF RESIDUALS TO POLYNOMIAL CLOCK MODELS 


The statistical characteristics of fractional frequency error and its 
integrated effect on range and Doppler observations have been discussed 
in detail. For range observations assume that the total random error 
is due to three sources, two of which are correlated noise processes. 
Then the total random range error is expressed as 

n(t) = Hg(t) + rij^(t) + C(t) (39) 

where q and are the correlated random range errors due to satel- 
lite an^ receiver random clock errors respectively. The quantity 4 
represents receiver white noise. The total integrated Doppler random 

error over the integration interval [t.,t„] is 

J * 

an(tj) = - n,(tj) t + t, (4o) 
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where is the white noise associated with the Doppler measurement 
procedure . 

Depending on the stability of the clock, the random range or Dop- 
pler error components, HgCi-) ^nd may appear quite systematic 

over fixed time intervals and may be represented by polynomial models 
of varying degree. For short time intervals the models for clock error 
were taken to be a bias and drift for range observations and a frequen- 
cy bias for Doppler observations. However, these models and even 
higher order polynomial models are not sufficient to entirely represent 
this correlated error. Thus knowledge of the statistical properties of 
the deviations of the error from such a model becomes important, as 
these residuals represent an unmodeled part of the observation equation 
after the inclusion of the polynomial model. 

Proceeding, equation (39) is expressed as follows 

where is an m'th degree polynomial chosen to model the correlat- 
ed random error Hg(t) and is an n'th degree polynomial model- 

ing the random process r)j^(t). The statistics of the range residuals 
r(t) may be developed from the covariance of the random clock errors. 
The second order statistics of the range residuals r(t) to a polynomial 
model are obtained as 

E[r(t)r''’(t)] = GE[R(t)R‘^(t)]G'^ (42) 


where 


G = [1 - A(A^A)‘^A^] (43) 

and A is the least squares design matrix for the polynomial model 

T 

selected. The E[R(t)R (t)] is the covariance matrix of the random 
clock error being modeled. This covariance is given by equations (32) 
and (33) . 

For integrated Doppler observations the statistics of the residuals to 
a given degree polynomial model are similarly obtained from equations 
(42) and (43), using the covariance matrix for integrated Doppler 
random error due to each system clock, equations (37) and (38). The 
equation may be written as 

E[Ar(t)Ar^(t)] = HE [AR(t)AR^(t) ]H^ (44) 
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where the matrix H is similar to the matrix G of equation (43) with 
changes due to the choice of the model adopted for clock-induced random 
Doppler errors 


T T 

H = [I - A'(A'A')'^A'] . (45) 
After the selection of the polynomial model, equation (40) has the form 

If the statistics of these residuals were ignored in an estimation 
problem, then the resulting parameter covariance matrix would be opti- 
mistic. An increase in the degrees of the polynomial clock models would 
offset this optimism to some extent since the level of unmodeled error 
would be decreased. However, if a rigorous estimation is to be per- 
formed, then these residual statistics must be included in the weight 
matrix to account for the unmodeled error r(t) or Ar(t) in a statistic- 
al rather than parametric fashion. The estimation algorithm should 
then produce a valid parameter covariance matrix regardless of the 
order of the polynomial models used provided numerical problems are not 
encountered and the parameters are independent and well observed. 

Finally, the theoretical standard errors for range residuals to a 
linear fit were determined using equation (42) for the rubidium and 
cesium clocks. The results are given in Figures 15 and 16. These 
figures graphically demonstrate that the statistics of the residuals to 
the clock modeling polynomial are not stationary. The variance of a 
residual depends on the order of the polynomial, the interval length 
and the location within the sample. However, the statistics of the 
residuals will be constant from interval to interval of the same length 
provided the sampling is performed equivalently and the same order 
polynomial is used. 
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Fig. 15-Standard error of satellite rubidium 
clock residuals based on a linear fit 



Fig. 16-Standard error of cesium clock 
residuals based on a linear fit 


576 


REFERENCES 


Allan, D. W. , J. H. Shoaf, and D. Halford, 1974. "Statistics of Time 
and Frequency Data Analysis," Time and Frequency; Theory and 
Fundamentals , National Bureau of Standards Monograph 140, Blair 
(Editor) . 

Anderle, R. J. , 1978. "Clock Performance as a Critical Parameter in 

Navigation Satellite Systems," Proceedings of the Tenth Annual 
Precise Time and Time Interval (PTTI) Applications and Planning 
Meeting , National Aeronautics and Space Administration Technical 
Memorandum 80250, Greenbelt, Maryland. 

Barnes, J. A. and S. Jarvis, 1971. "Efficient Numerical and Analog 
Modeling of Flicker Noise Processes," National Bureau of Standards 
Technical Note 604, Boulder, Colorado. 

Barnes, J. A., A. R. Chi, L. S. Cutler, D. J. Healey, D. B. Leeson, 
T. E. McGunigal, J. A. Mullen, Jr., W. L. Smith, R. L. Sydnor, 
R. F. C. Vessot, and G. M. R. Winkler, 1971. "Characterization of 
Frequency Stability," IEEE Transactions on Instrument and Measure- 
ments , Vol. 20, No. 2. 

Bartholomew, C. A., 1978. "Satellite Frequency Standards," Navigation , 
Vol. 25, No. 2. 

Blair, B. E. (Editor), 1974. Time and Frequency: Theory and Fundamen- 

tals , U. S. Department of Commerce/National Bureau of Standards 
Monograph 140, U. S. Government Printing Office, Washington, D.C. 

Fell, P. J. and B. R. Hermann, 1979. "An Empirical Evaluation of the 
Effect of Oscillator Errors on Dynamic Point Positioning Based on 
the NAVSTAR GPS System," Proceedings of the Second International 
Geodetic Symposium on Satellite Doppler Positioning , Applied Re- 
search Laboratories, The University of Texas at Austin. 

Fell, P. J., 1980. "Geodetic Positioning Using a Global Positioning 

System of Satellites," Department of Geodetic Science Report No. 
299, The Ohio State University, Columbus. 

Gelb, A., (Editor), 1974. Applied Optimal Estimation , MIT Press, 
Cambridge, Massachusetts. 

Hellwig, H. W. , 1975. "Atomic Frequency Standards: A Survey," Pro- 

ceedings of the IEEE , Vol. 63, No. 2. 


577 


Hellwig, H. W. , 1977. "Frequency Standards and Clocks: A Tutorial 

Introduction," National Bureau of Standards Technical Note 616, 
U. S. Government Printing Office, Washington, D.C. 

Meditch, J. S. , 1975. "Clock Error Models for Simulation and Estima- 
tion," The Aerospace Corporation, Report No. TOR-0076 (647A-01)-2, 
El Segundo, California. 


578 


QUESTIONS AND ANSWERS 


DR. KARTASCHOFF: 

I have just one question that you have been remodeling the flicker 
noise level with these five processes, and I was asking myself now 
when I hear it what do you think about if one just could use direct- 
ly the time interval error estimation, as I had shown just before, 
also for estimating the range error? 

Furthermore, one could try to estimate the uncertainty of that 
estimation by using the uncertainty of the Allan variance using the 
theory of Audoin-Lesage that limited the sample that we always, for 
a given time, t, on the flicker level, and we always have an un- 
certainty that is given as the number of samples. For the last 
point you measured, you have only two samples; so you have 2000 per- 
cent error as the uncertainty. 

I think it would be an interesting exercise to repeat the cal- 
culation using these estimations and using your process. Very 
probably both will give very similar results and both can be used. 
That would be interesting, I think. 

MR. FELL: 

I think you are right. This is just one way that you could do this 
approximate method, those Allan variances. They are specified for 
a clock and are only an approximation of actual performance, but for 
■ a long time people have ignored this type of residual error which is 
left in the estimation problem. And I think that because we are now 
trying to get down to such small errors, namely baseline errors of 
less than 10 centimeters, that we are going to have to take a second 
look at our modeling and make sure that it is sufficient for the 
problems that we are addressing. 

Otherwise the parameter of statistics which we get out of the 
estimation algorithm are going to be too optimistic. 

DR. VICTOR REINHARDT, NASA/Goddard 

You also, if you do your conversions to your range statistics, or 
range rate statistics, up front you will find that your range rate 
estimator is one of the weighted Allan variances, and you can just 
use the tables that are published from converting to the zero dead 
time Allan variance to the Allan variance with dead times to get 
your range statistics. That all combined with some factor like C 

or-C . 
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And you can do this not even with a hand calculator on the 
back of an envelope, do the same thing just by looking up the NBS 
publications on the various weighings for the various models since 
all of the frequency standards that we use breakup into well de- 
fined regimes we have, you know, a well defined parallel. 
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